d <- read.csv("rbalpha.csv")[, -1]
d[d == 8] <- NA
d[d == 7] <- 0
d <- as.matrix(d)
d1 <- d[1:9, ]
d2 <- d[10:18, ]
d3 <- d[19:27, ]
d4 <- d[28:36, ]
d5 <- d[37:45, ]
d1[upper.tri(d1)] <- t(d1)[upper.tri(d1)]
d2[upper.tri(d2)] <- t(d2)[upper.tri(d2)]
d3[upper.tri(d3)] <- t(d3)[upper.tri(d3)]
d4[upper.tri(d4)] <- t(d4)[upper.tri(d4)]
d5[upper.tri(d5)] <- t(d5)[upper.tri(d5)]

balphas <- simplify2array(list(d1, d2, d3, d4, d5))

balpha <- function(q1, q0, verbose = FALSE) {
	if(q1 < 4 | q0 < 4) stop("number of clusters too small")
	ba <- balphas[q1-3, q0-3, ]
	names(ba) <- c("10%", "5%", "2.5%", "1%", "0.5%")
	if (verbose == TRUE) {
		ba[is.na(ba)] <- "not available"
		ba[ba == 0] <- "use second largest order statistic"
		cat(" ", names(ba)[1], ": ", ba[1] , "\n", sep = "")
		cat("  ", names(ba)[2], ": ", ba[2] , "\n", sep = "")
		cat(names(ba)[3], ": ", ba[3] , "\n", sep = "")
		cat("  ", names(ba)[4], ": ", ba[4] , "\n", sep = "")
		cat(names(ba)[5], ": ", ba[5] , "\n", sep = "")
	} else {
		ba
	}
}

MakeIndex <- function(q1, q0) {
	count.stats <- 1:(q1  + q0)
	index.treat <- combn(count.stats, q1)
	index.cont <- sapply(1:ncol(index.treat), function(j) count.stats[! count.stats %in% index.treat[, j]])	
	index <- cbind(t(index.treat), t(index.cont))
	index
}

MakeRandIndex <- function(q1, q0, n.greps = 9999) {
	q <- q1 + q0
	index <- t(cbind(1:q, replicate(n.greps, sample(1:q, q))))
	index
}

OrdTest <- function(stats, k, index, q1, q0) {
	q <- q1 + q0
	M <- nrow(index)
	r.stats <- matrix(stats[index[1:M, ]], M, q)
	perm.stats <- rowMeans(r.stats[, 1:q1]) - rowMeans(r.stats[, (q1+1):q])
	orig.stat <- perm.stats[1]
	perm.stats <- sort(perm.stats)
	as.numeric(orig.stat > perm.stats[k])
}

IMTest <- function(stats, alpha = .05, q1, q0) {
	q <- q1 + q0
	t.stat <- (mean(stats[1:q1]) - mean(stats[(q1+1):q]))/sqrt(var(stats[1:q1])/q1 + var(stats[(q1+1):q])/q0)
	as.numeric(t.stat > qt(1-alpha, min(q1, q0) - 1))
}


sims <- 10^4
q1 <- 8
q0 <- 8
idx <- MakeIndex(q1,q0)
#idx <- MakeRandIndex(q1,q0)
sdmat <- matrix( 
	c(rep(1, q1 + q0),
	  rep(1, q1 + q0-1), rep(10^2, 1),
	  rep(1, q1-1), rep(10^2, 1), rep(1, q0-1), rep(10^2, 1),
	  rep(1, q1/2), rep(1, q1/2), rep(3, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(1, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(3, q0/2), rep(9, q0/2)
	  ), 6, q1 + q0, TRUE
	  )

muvec <- c(2, 40, 40, 5, 5, 10)

rs.ap <- matrix(NA, 6, 51)
rs.im <- rs.ap

tst <- function(x) c(OrdTest(x, ceiling(round(((1-balpha(q1, q0)[2])*dim(idx)[1]), 4)), idx, q1, q0), IMTest(x, .05, q1, q0))

for(i in 1:6) {
	mus <- seq(0, muvec[i], length.out = 51)
	n.mus <- length(mus)
	sds <- sdmat[i, ]
	res <- matrix(NA, 2, n.mus)
	for(s in 1:n.mus) {
		mu <- c(rep(mus[s], q1), rep(0,q0))
		set.seed(48109)
		res[, s] <- rowMeans(replicate(sims, tst(rnorm(q1+q0)*sds + mu)))
	}
	rs.ap[i, ] <- res[1, ]
	rs.im[i, ] <- res[2, ]
	print(rs.ap)
	print(rs.im)
}

res_norm_8_8 <- rbind(rs.ap, rs.im)

write.table(rbind(rs.ap, rs.im), "rperm_norm_loc_rev_8_8.txt")


q1 <- 8
q0 <- 12
#idx <- MakeIndex(q1,q0)
idx <- MakeRandIndex(q1,q0)
sdmat <- matrix( 
	c(rep(1, q1 + q0),
	  rep(1, q1 + q0-1), rep(10^2, 1),
	  rep(1, q1-1), rep(10^2, 1), rep(1, q0-1), rep(10^2, 1),
	  rep(1, q1/2), rep(1, q1/2), rep(3, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(1, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(3, q0/2), rep(9, q0/2)
	  ), 6, q1 + q0, TRUE
	  )

muvec <- c(2, 40, 40, 5, 5, 10)

rs.ap <- matrix(NA, 6, 51)
rs.im <- rs.ap

tst <- function(x) c(OrdTest(x, ceiling(round(((1-balpha(q1, q0)[2])*dim(idx)[1]), 4)), idx, q1, q0), IMTest(x, .05, q1, q0))

for(i in 1:6) {
	mus <- seq(0, muvec[i], length.out = 51)
	n.mus <- length(mus)
	sds <- sdmat[i, ]
	res <- matrix(NA, 2, n.mus)
	for(s in 1:n.mus) {
		mu <- c(rep(mus[s], q1), rep(0,q0))
		set.seed(48109)
		res[, s] <- rowMeans(replicate(sims, tst(rnorm(q1+q0)*sds + mu)))
	}
	rs.ap[i, ] <- res[1, ]
	rs.im[i, ] <- res[2, ]
	print(rs.ap)
	print(rs.im)
}

res_norm_8_12 <- rbind(rs.ap, rs.im)

write.table(rbind(rs.ap, rs.im), "rperm_norm_loc_rev_8_16.txt")


q1 <- 16
q0 <- 16
#idx <- MakeIndex(q1,q0)
idx <- MakeRandIndex(q1,q0)
sdmat <- matrix( 
	c(rep(1, q1 + q0),
	  rep(1, q1 + q0-1), rep(10^2, 1),
	  rep(1, q1-1), rep(10^2, 1), rep(1, q0-1), rep(10^2, 1),
	  rep(1, q1/2), rep(1, q1/2), rep(3, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(1, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(3, q0/2), rep(9, q0/2)
	  ), 6, q1 + q0, TRUE
	  )

muvec <- c(2, 40, 40, 5, 5, 10)

rs.ap <- matrix(NA, 6, 51)
rs.im <- rs.ap

tst <- function(x) c(OrdTest(x, ceiling(round(((1-.0491)*dim(idx)[1]), 4)), idx, q1, q0), IMTest(x, .05, q1, q0))

for(i in 1:6) {
	mus <- seq(0, muvec[i], length.out = 51)
	n.mus <- length(mus)
	sds <- sdmat[i, ]
	res <- matrix(NA, 2, n.mus)
	for(s in 1:n.mus) {
		mu <- c(rep(mus[s], q1), rep(0,q0))
		set.seed(48109)
		res[, s] <- rowMeans(replicate(sims, tst(rnorm(q1+q0)*sds + mu)))
	}
	rs.ap[i, ] <- res[1, ]
	rs.im[i, ] <- res[2, ]
	print(rs.ap)
	print(rs.im)
}

res_norm_16_16 <- rbind(rs.ap, rs.im)

write.table(rbind(rs.ap, rs.im), "rperm_norm_loc_rev_16_16.txt")



sims <- 10^4
q1 <- 8
q0 <- 8
idx <- MakeIndex(q1,q0)
#idx <- MakeRandIndex(q1,q0)
sdmat <- matrix( 
	c(rep(1, q1 + q0),
	  rep(1, q1 + q0-1), rep(10^2, 1),
	  rep(1, q1-1), rep(10^2, 1), rep(1, q0-1), rep(10^2, 1),
	  rep(1, q1/2), rep(1, q1/2), rep(3, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(1, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(3, q0/2), rep(9, q0/2)
	  ), 6, q1 + q0, TRUE
	  )

muvec <- c(2, 40, 40, 5, 5, 10)

rs.ap <- matrix(NA, 6, 51)
rs.im <- rs.ap

tst <- function(x) c(OrdTest(x, ceiling(round(((1-balpha(q1, q0)[2])*dim(idx)[1]), 4)), idx, q1, q0), IMTest(x, .05, q1, q0))

for(i in 1:6) {
	mus <- seq(0, muvec[i], length.out = 51)
	n.mus <- length(mus)
	sds <- sdmat[i, ]
	res <- matrix(NA, 2, n.mus)
	for(s in 1:n.mus) {
		mu <- c(rep(mus[s], q1), rep(0,q0))
		set.seed(48109)
		res[, s] <- rowMeans(replicate(sims, tst(rt(q1+q0, 1)*sds + mu)))
	}
	rs.ap[i, ] <- res[1, ]
	rs.im[i, ] <- res[2, ]
	print(rs.ap)
	print(rs.im)
}

res_cauchy_8_8 <- rbind(rs.ap, rs.im)

write.table(rbind(rs.ap, rs.im), "rperm_cauchy_loc_rev_8_8.txt")


q1 <- 8
q0 <- 12
#idx <- MakeIndex(q1,q0)
idx <- MakeRandIndex(q1,q0)
sdmat <- matrix( 
	c(rep(1, q1 + q0),
	  rep(1, q1 + q0-1), rep(10^2, 1),
	  rep(1, q1-1), rep(10^2, 1), rep(1, q0-1), rep(10^2, 1),
	  rep(1, q1/2), rep(1, q1/2), rep(3, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(1, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(3, q0/2), rep(9, q0/2)
	  ), 6, q1 + q0, TRUE
	  )

muvec <- c(2, 40, 40, 5, 5, 10)

rs.ap <- matrix(NA, 6, 51)
rs.im <- rs.ap

tst <- function(x) c(OrdTest(x, ceiling(round(((1-balpha(q1, q0)[2])*dim(idx)[1]), 4)), idx, q1, q0), IMTest(x, .05, q1, q0))

for(i in 1:6) {
	mus <- seq(0, muvec[i], length.out = 51)
	n.mus <- length(mus)
	sds <- sdmat[i, ]
	res <- matrix(NA, 2, n.mus)
	for(s in 1:n.mus) {
		mu <- c(rep(mus[s], q1), rep(0,q0))
		set.seed(48109)
		res[, s] <- rowMeans(replicate(sims, tst(rt(q1+q0, 1)*sds + mu)))
	}
	rs.ap[i, ] <- res[1, ]
	rs.im[i, ] <- res[2, ]
	print(rs.ap)
	print(rs.im)
}

res_cauchy_8_12 <- rbind(rs.ap, rs.im)

write.table(rbind(rs.ap, rs.im), "rperm_cauchy_loc_rev_8_16.txt")


q1 <- 16
q0 <- 16
#idx <- MakeIndex(q1,q0)
idx <- MakeRandIndex(q1,q0)
sdmat <- matrix( 
	c(rep(1, q1 + q0),
	  rep(1, q1 + q0-1), rep(10^2, 1),
	  rep(1, q1-1), rep(10^2, 1), rep(1, q0-1), rep(10^2, 1),
	  rep(1, q1/2), rep(1, q1/2), rep(3, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(1, q0/2), rep(3, q0/2),
	  rep(1, q1/2), rep(3, q1/2), rep(3, q0/2), rep(9, q0/2)
	  ), 6, q1 + q0, TRUE
	  )

muvec <- c(2, 40, 40, 5, 5, 10)

rs.ap <- matrix(NA, 6, 51)
rs.im <- rs.ap

tst <- function(x) c(OrdTest(x, ceiling(round(((1-.0491)*dim(idx)[1]), 4)), idx, q1, q0), IMTest(x, .05, q1, q0))

for(i in 1:6) {
	mus <- seq(0, muvec[i], length.out = 51)
	n.mus <- length(mus)
	sds <- sdmat[i, ]
	res <- matrix(NA, 2, n.mus)
	for(s in 1:n.mus) {
		mu <- c(rep(mus[s], q1), rep(0,q0))
		set.seed(48109)
		res[, s] <- rowMeans(replicate(sims, tst(rt(q1+q0, 1)*sds + mu)))
	}
	rs.ap[i, ] <- res[1, ]
	rs.im[i, ] <- res[2, ]
	print(rs.ap)
	print(rs.im)
}

res_cauchy_16_16 <- rbind(rs.ap, rs.im)

write.table(rbind(rs.ap, rs.im), "rperm_cauchy_loc_rev_16_16.txt")




pdf("rperm_norm_loc_rev.pdf", height = 10, width = 6)
par(mfrow = c(6, 3), mar=c(2, 4, 2, 1))

plot(seq(0, muvec[1], length.out = 51), res_norm_8_8[7, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(a)",side=2,las=1,line=2)
mtext(expression(paste(italic(q)[1], "= 8, ", italic(q)[0], "= 8")),side=3,las=1,line=.5)
lines(seq(0, muvec[1], length.out = 51), res_norm_8_8[1, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")
legend("topleft", legend = c("Adj. perm.", "IM", "5% lvl."), 
								lty =c(1, 1, 2), lwd = c(rep(1.5, 3)), 
								col = c("black", "grey", "black"),
								pt.cex = 1, cex = .9, bty = "n")

plot(seq(0, muvec[1], length.out = 51), res_norm_8_12[7, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext(expression(paste(italic(q)[1], "= 8, ", italic(q)[0], "= 16")),side=3,las=1,line=.5)
lines(seq(0, muvec[1], length.out = 51), res_norm_8_12[1, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[1], length.out = 51), res_norm_16_16[7, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext(expression(paste(italic(q)[1], "= 16, ", italic(q)[0], "= 16")),side=3,las=1,line=.5)
lines(seq(0, muvec[1], length.out = 51), res_norm_16_16[1, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[2], length.out = 51), res_norm_8_8[8, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(b)",side=2,las=1,line=2)
lines(seq(0, muvec[2], length.out = 51), res_norm_8_8[2, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[2], length.out = 51), res_norm_8_12[8, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[2], length.out = 51), res_norm_8_12[2, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[2], length.out = 51), res_norm_16_16[8, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[2], length.out = 51), res_norm_16_16[2, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[3], length.out = 51), res_norm_8_8[9, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(c)",side=2,las=1,line=2)
lines(seq(0, muvec[3], length.out = 51), res_norm_8_8[3, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[3], length.out = 51), res_norm_8_12[9, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[3], length.out = 51), res_norm_8_12[3, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[3], length.out = 51), res_norm_16_16[9, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[3], length.out = 51), res_norm_16_16[3, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[4], length.out = 51), res_norm_8_8[10, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(d)",side=2,las=1,line=2)
lines(seq(0, muvec[4], length.out = 51), res_norm_8_8[4, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[4], length.out = 51), res_norm_8_12[10, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[4], length.out = 51), res_norm_8_12[4, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[4], length.out = 51), res_norm_16_16[10, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[4], length.out = 51), res_norm_16_16[4, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[5], length.out = 51), res_norm_8_8[11, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(e)",side=2,las=1,line=2)
lines(seq(0, muvec[5], length.out = 51), res_norm_8_8[5, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[5], length.out = 51), res_norm_8_12[11, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[5], length.out = 51), res_norm_8_12[5, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[5], length.out = 51), res_norm_16_16[11, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[5], length.out = 51), res_norm_16_16[5, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[6], length.out = 51), res_norm_8_8[12, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(f)",side=2,las=1,line=2)
lines(seq(0, muvec[6], length.out = 51), res_norm_8_8[6, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[6], length.out = 51), res_norm_8_12[12, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[6], length.out = 51), res_norm_8_12[6, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[6], length.out = 51), res_norm_16_16[12, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[6], length.out = 51), res_norm_16_16[6, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

dev.off()




pdf("rperm_cauchy_loc_rev.pdf", height = 10, width = 6)
par(mfrow = c(6, 3), mar=c(2, 4, 2, 1))

plot(seq(0, muvec[1], length.out = 51), res_cauchy_8_8[7, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(a)",side=2,las=1,line=2)
mtext(expression(paste(italic(q)[1], "= 8, ", italic(q)[0], "= 8")),side=3,las=1,line=.5)
lines(seq(0, muvec[1], length.out = 51), res_cauchy_8_8[1, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")
legend("topleft", legend = c("Adj. perm.", "IM", "5% lvl."), 
								lty =c(1, 1, 2), lwd = c(rep(1.5, 3)), 
								col = c("black", "grey", "black"),
								pt.cex = 1, cex = .9, bty = "n")

plot(seq(0, muvec[1], length.out = 51), res_cauchy_8_12[7, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext(expression(paste(italic(q)[1], "= 8, ", italic(q)[0], "= 16")),side=3,las=1,line=.5)
lines(seq(0, muvec[1], length.out = 51), res_cauchy_8_12[1, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[1], length.out = 51), res_cauchy_16_16[7, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext(expression(paste(italic(q)[1], "= 16, ", italic(q)[0], "= 16")),side=3,las=1,line=.5)
lines(seq(0, muvec[1], length.out = 51), res_cauchy_16_16[1, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[2], length.out = 51), res_cauchy_8_8[8, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(b)",side=2,las=1,line=2)
lines(seq(0, muvec[2], length.out = 51), res_cauchy_8_8[2, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[2], length.out = 51), res_cauchy_8_12[8, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[2], length.out = 51), res_cauchy_8_12[2, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[2], length.out = 51), res_cauchy_16_16[8, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[2], length.out = 51), res_cauchy_16_16[2, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[3], length.out = 51), res_cauchy_8_8[9, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(c)",side=2,las=1,line=2)
lines(seq(0, muvec[3], length.out = 51), res_cauchy_8_8[3, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[3], length.out = 51), res_cauchy_8_12[9, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[3], length.out = 51), res_cauchy_8_12[3, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[3], length.out = 51), res_cauchy_16_16[9, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[3], length.out = 51), res_cauchy_16_16[3, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[4], length.out = 51), res_cauchy_8_8[10, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(d)",side=2,las=1,line=2)
lines(seq(0, muvec[4], length.out = 51), res_cauchy_8_8[4, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[4], length.out = 51), res_cauchy_8_12[10, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[4], length.out = 51), res_cauchy_8_12[4, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[4], length.out = 51), res_cauchy_16_16[10, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[4], length.out = 51), res_cauchy_16_16[4, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[5], length.out = 51), res_cauchy_8_8[11, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(e)",side=2,las=1,line=2)
lines(seq(0, muvec[5], length.out = 51), res_cauchy_8_8[5, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[5], length.out = 51), res_cauchy_8_12[11, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[5], length.out = 51), res_cauchy_8_12[5, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[5], length.out = 51), res_cauchy_16_16[11, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[5], length.out = 51), res_cauchy_16_16[5, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[6], length.out = 51), res_cauchy_8_8[12, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
mtext("(f)",side=2,las=1,line=2)
lines(seq(0, muvec[6], length.out = 51), res_cauchy_8_8[6, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[6], length.out = 51), res_cauchy_8_12[12, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[6], length.out = 51), res_cauchy_8_12[6, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

plot(seq(0, muvec[6], length.out = 51), res_cauchy_16_16[12, ], ylim = c(0,1), type = "l", lwd = 1.5, col = "grey", ylab = "", xlab = "")
lines(seq(0, muvec[6], length.out = 51), res_cauchy_16_16[6, ], lwd = 1.5)
abline(.05, 0, lwd = 1.5, lty = "dashed")

dev.off()




